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1.  Introduction 


Traumatic  brain  injury  (TBI)  is  a  serious  concern  for  the  military  and  the  general 
civilian  population.  Blast-related  TBI  has  been  prevalent  in  recent  military  conflicts 
(Gupta  and  Przekwas  2013).  During  Operation  Iraqi  Freedom,  a  study  of  casualties 
requiring  level  V  care  at  Walter  Reed  Army  Medical  Center  reported  that  29%  of 
those  that  were  screened  had  a  TBI.  Blasts  and  explosions  were  the  most  common 
causes,  accounting  for  78%  of  those  found  to  have  a  TBI  (Traumatic  Brain  Injury 
Task  Force  2007).  Understanding  TBI  also  has  significant  relevance  for  the  civilian 
population.  Approximately  1.4  million  people  within  the  United  States  sustain  a 
TBI  each  year.  Of  that  number  50,000  die,  235,000  are  hospitalized,  and  1 . 1  million 
are  evaluated,  treated  and  released  from  emergency  departments  (Langlois  et  al. 
2006).  When  one  also  considers  concussions  (often  called  mild  TBI)  it  is  possible 
the  largest  proportion  of  patients  are  not  even  seen  in  an  emergency  department. 
TBI  is  even  more  concerning  due  to  its  residual  effects.  It  is  estimated  that  at  least 
5.3  million  Americans,  almost  2%  of  the  population,  have  current  long-term  or 
lifelong  disabilities  as  a  result  of  TBI  (Defense  and  Veterans  Brain  Injury  Center 
2014). 

TBI  associated  with  closed  head  injuries,  also  referred  to  as  nonpenetrating  head 
injuries,  can  be  caused  by  blast,  blunt-force  impact,  or  sudden  acceleration.  In  cases 
such  as  these,  diffuse  axonal  injury  is  one  particular  injury  mechanism  that  has  been 
cited  as  a  signature  injury  of  TBI  neural  damage  (Taber  et  al.  2006;  Gupta  and 
Przekwas  2013).  Defonnation  of  the  brain  tissue  can  induce  misalignment  in  the 
cytoskeletal  network  or  axolemma  permeability,  inducing  a  cascade  of  subcellular 
events  culminating  in  the  severance  of  the  axon  (Christman  et  al.  1994;  Smith  et  al. 
2003).  It  is  these  axon  fiber  bundles  that  make  up  the  structural  network  that  allows 
neurons  to  communicate  with  one  another.  Injury  to  the  axons  leads  to  degraded 
structural  connectivity,  which  may  be  responsible  for  the  cognitive  deficits  that  are 
characteristic  of  mild,  moderate,  and  severe  cases  of  TBI  (Vettel  et  al.  2010). 
Concussion,  or  mild  TBI,  is  thought  to  be  a  less  severe  type  of  diffuse  axonal  injury 
where  axons  are  damaged  to  a  minor  extent  from  stretching.  Postmortem  studies  of 
brains  with  concussions  have  found  axonal  damage;  however,  because  of  other 
factors,  such  as  restricted  blood  flow,  it  is  not  possible  to  isolate  the  cause  of  this 
damage  and  solely  link  it  to  concussions. 

While  substantial  work  has  already  been  performed  by  the  Soldier  Protection 
Sciences  Branch  to  investigate  axonal  injury,  another  aspect  of  TBI,  hemorrhage 
and  other  forms  of  vascular  injury,  has  not  been  represented  within  the  numerical 
models.  In  2006  a  study  of  head  injury  sustained  by  Soldiers  in  Operation  Iraqi 
Freedom  was  completed  by  Dr  Rocco  Armonda.  Of  the  patients  in  the  study,  82% 
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suffered  blast  injury  and  14%  were  injured  by  a  gunshot  wound.  Of  the  entire  group, 
injury  to  the  cerebral  vasculature  was  seen  in  a  significant  number  of  patients.  The 
data  indicate  that  47%  had  vasospasm,  35%  had  a  pseudoaneurysm,  12%  had  a 
subarachnoid  hemorrhage,  3.5%  had  an  epidural  hematoma,  16%  had  a  subdural 
hematoma,  11%  had  an  intraventricular  hemorrhage,  and  14%  had  mixed 
hemorrhages  (Armonda  et  al.  2006).  In  order  for  our  numerical  models  to  better 
predict  the  extent  of  injury  to  the  Solider,  a  method  to  include  the  cerebral 
vasculature  must  be  implemented.  This  report  details  one  such  approach,  which 
builds  upon  existing  simulation  work  within  the  Soldier  Protection  Sciences 
Branch,  and  through  a  postprocessing  algorithm,  determines  the  strains  within  the 
cerebral  vascular  network. 

2.  Procedure 


The  following  approach  to  modeling  injury  to  the  cerebral  vasculature  takes 
advantage  of  an  existing  head  and  brain  model  developed  within  the  Soldier 
Protection  Sciences  Branch  (Kraft  et  al.  2010;  Kraft  and  Dagro  2011;  Kraft  et  al. 
2012;  Dagro  et  al.  2013;  McKee  et  al.  2013).  The  algorithm  will  use  the  results 
from  those  simulations  to  calculate  the  deformation  and  strain  within  a  cerebral 
vascular  network.  The  intention  is  to  leave  the  existing  head  and  brain  mesh 
untouched  and  create  an  overlapping  beam  mesh  of  the  cerebral  vasculature 
(Fig.  1).  The  proposed  algorithm  would  serve  as  a  link  between  the  2  meshes  such 
that  deformation  and  strain  could  be  mapped  from  the  full  finite  element  simulation 
of  the  head  and  brain  onto  vascular  elements.  Thus,  the  vascular  model  would  not 
constitute  a  full  finite  element  simulation  on  its  own  but  rather  a  postprocessing  of 
the  results  from  an  already  existing  finite  element  model.  The  deformation  data 
mapped  onto  the  vascular  network  can  then  be  fed  into  a  damage  model  selected  by 
the  user  to  identify  the  locations  at  greatest  risk  for  injury. 
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Fig.  1  A  3-dimensional  representation  of  the  cerebral  vasculature  showing  larger  arteries. 
Smaller  capillaries  are  not  pictured  but  can  also  be  included. 

This  approach  would  require  no  alternations  to  the  existing  simulation.  However, 
once  the  algorithm  has  detennined  the  linkages  between  the  2  meshes,  it  is  possible 
to  use  the  information  regarding  the  location  and  orientation  of  the  vascular 
segments  to  generate  an  updated  material  model  for  the  elements  within  the  brain. 
The  presence  of  the  vasculature  running  through  the  brain  provides  an  increased 
stiffness  to  its  mechanical  response.  In  most  simulations  this  is  not  directly 
modeled,  instead  the  stiffness  of  the  brain  tissue  is  increased  to  account  for  the 
vascular  presence.  Unfortunately,  such  an  approach  fails  to  capture  the  directional 
nature  of  the  vascular  network.  Therefore,  by  using  the  location  of  the  vasculature 
and  its  orientation  a  more  accurate  anisotropic  material  model  can  be  employed  in 
the  initial  head  and  brain  simulation.  As  before,  this  would  still  be  separate  from 
the  vascular  model,  only  with  new  material  properties  defined  for  the  elements 
containing  vascular  segments. 

Thus,  this  project  will  provide  both  a  pre-  and  a  postprocessing  algorithm  built 
around  the  linking  of  a  solid  element  brain  mesh  with  a  beam  network  mesh  for  the 
cerebral  vasculature.  The  following  subsections  will  begin  by  detailing  the  linking 
of  the  2  independent  meshes,  the  common  core  of  the  algorithm,  before  moving  on 
to  the  pre-  and  postprocessing  portions  of  the  code. 

2.1  Reading  Mesh  Information  and  Output  Data 

To  determine  the  links  between  the  solid  element  brain  mesh  and  the  vascular 
network  mesh,  and  manipulate  the  deformation  data  produced  by  the  finite  element 
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simulation,  the  proposed  algorithm  must  be  able  to  interface  with  the  files  generated 
by  the  meshing  programs  and  the  results  produced  by  finite  element  codes  such  as 
SIERRA.  Reading  keyword  format  input  files  such  as  those  used  in  commercial 
codes  like  LS-DYNA  is  relatively  straightforward  and  can  be  done  using  the 
standard  input  and  output  libraries  for  C/C++.  However,  other  finite  element  codes, 
including  SIERRA,  use  the  EXODUS  II  file  format  for  both  mesh  data  and  output 
data.  The  EXODUS  II  file  fonnat  was  developed  at  Sandia  National  Laboratories 
specifically  to  store  and  retrieve  data  for  finite  element  analyses.  To  interface  with 
the  files  in  the  EXODUS  II  format,  special  libraries  and  commands  are  required 
that  are  not  included  in  standard  C/C++  packages. 

The  files  necessary  to  generate  the  EXODUS  II  libraries  are  available  and  can  be 
retrieved  from  sites  like  SourceForge.net.  The  user  manual  is  included  in 
Appendix  A.  A  sample  input  file  is  provided  in  Appendix  B  and  alternate  fonnat 
for  postprocessing  data  is  provided  in  Appendix  C.  A  list  of  useful  EXODUS  II 
commands  is  included  in  Appendix  D. 

2.2  Relating  Vessel  Locations  to  Brain  Mesh 

Since  the  initial  locations  of  the  nodes  and  elements  associated  with  the  brain  mesh 
are  independent  of  the  location  of  the  nodes  and  elements  associated  with  the 
vascular  mesh,  an  algorithm  is  needed  to  relate  the  2.  It  is  necessary  to  know  which 
solid  (brain)  element  contains  a  given  beam  (vascular)  element  to  calculate  the 
strains  in  the  vascular  mesh,  and  use  the  orientation  of  the  vasculature  to  infonn  the 
anisotropic  response  of  the  surrounding  brain  tissue.  However,  because  most  of  the 
calculations  are  done  in  an  element  reference  domain,  as  opposed  to  the  physical 
(x,  y,  z)  domain,  it  is  also  necessary  to  relate  a  given  point  on  a  beam  element 
(x,  y,  z)  to  the  reference  coordinates  (<f,  rj,  Q  for  the  surrounding  solid  element  as 
shown  in  Fig.  2. 
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Fig.  2  Mapping  from  the  current  domain  (x,  y,  z-space)  to  the  element  domain  (f,  77, 
space) 

Mapping  a  point  from  the  element  domain  to  the  physical  domain  can  be  done  with 
the  element  shape  functions,  Na,  and  the  location  of  the  nodes  in  the  physical 
domain,  X_a,  using  the  following  equation: 


NEN 


where  NEN  is  the  number  of  nodes  and  shapes  functions  associated  with  that 
element.  For  a  typical  simulation  of  the  brain  and  head,  the  solid  elements  will  be 
either  trilinear  hexahedral  elements  or  linear  tetrahedral  elements.  The  shape 
functions  for  a  trilinear  hexahedral  element  are 

Natf.V.O  =  5(1  +  faO(l  +  VaVX  1  +  (aO,  (2) 

where  the  parameters  for  the  shape  functions  are  given  in  Table  1. 


Table  1  Hexahedral  shape  function  parameters 


a 

V  a 

1 

-1 

-1 

1 

2 

1 

-1 

1 

3 

1 

1 

1 

4 

-1 

1 

1 

5 

-1 

-1 

-1 

6 

1 

-1 

-1 

7 

1 

1 

-1 

8 

-1 

1 

-1 
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The  shape  functions  for  the  linear  tetrahedral  elements  are  as  follows: 


(3) 

n2(£>v>0  =  v , 

(4) 

N&.V.O  =  <■ ,  and 

(5) 

(6) 

While  the  shape  functions  are  a  convenient  way  to  map  from  (<f ,  tj,  £)  to  (x,  y,  z), 
what  is  required  here  is  the  reverse  mapping  (x,  y,  z)  to  (<f,  rj,  <').  While  this  can  be 
derived  for  a  linear  tetrahedral  element,  a  solution  cannot  be  simply  written  for  a 
general  trilinear  hexahedral  element.  Instead  we  must  employ  a  solution  technique 
such  as  Newton’s  method  to  solve  the  nonlinear  system  of  equations.  Because  this 
approach  is  computationally  intensive,  the  search  algorithm  is  broken  up  into  an 
initial  coarse  search,  which  is  inexpensive  to  perfonn,  followed  by  a  fine  search 
where  Newton’s  method  is  employed.  For  large  meshes,  the  reduction  in 
computational  time  can  be  significant. 

Prior  to  initiating  the  coarse  search,  the  location  of  centroid,  X_c ,  for  each  solid 
element  is  calculated 


l 

NEN 


NEN 


I 


xa 


(7) 


Following  that,  the  radius  of  the  smallest  sphere  containing  all  points  within  a  given 
element  and  centered  at  the  element’s  centroid  is  calculated  (Fig.  3). 


Fig.  3  The  centroid  of  the  element  and  sphere  just  large  enough  to  contain  all  of  the  points 
within  the  element 

To  perform  the  coarse  search,  the  distance  from  the  centroid  of  each  solid  element 
to  a  given  point  in  the  vascular  mesh  is  calculated.  If  that  distance  is  less  than  or 
equal  to  the  previously  calculated  radius  associated  with  that  solid  element,  then 
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the  element  potentially  contains  that  point  and  a  fine  level  search  is  required.  If  the 
distance  to  the  point  is  greater  than  the  radius,  then  the  point  cannot  fall  within  the 
element,  and  the  fine  level  search  is  unnecessary. 


The  fine  level  search  will  determine  whether  the  point  truly  falls  within  the  element, 
like  point  p2  in  Fig.  4,  or  if  it  is  still  outside,  like  point  pi.  To  do  this  we  must  solve 
the  following  system  of  equations  for  77, 

NEN 

Y  NM.r],Ola  ~X  =  0 

,  (8) 

where  X_a  are  the  nodal  positions  for  the  given  element  and  X_  is  the  location  of  the 
vascular  point  in  the  physical  domain,  both  known  quantities.  To  begin  we  rewrite 
this  as  a  function,  f  —  [Safi  Na  (■fj  ^aJ  —X_—0,  where  the  coordinates 
rj,  (  have  been  rewritten  as  the  vector  .  Without  knowing  the  value  of  that 
satisfies  the  equation,  we  can  make  an  initial  guess,  °,  and  Taylor  expand  about  it 


0  = 


(9) 


Fig.  4  Points  contained  within  the  sphere  can  also  be  contained  by  the  element,  although  it 
is  not  necessarily  so 


Dropping  higher  order  terms,  we  can  solve  for  This  new  solution  can  be 

substituted  for  the  previous  guess  and  process  repeated,  giving  us  the  following 
equation: 


n+l 


(10) 


The  gradient  of  /  can  be  calculated  from  the  shape  functions  and  the  nodal 


positions 
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a_M 

dl 


NEN 


2> 


w-(f) 

dl 


(11) 


This  is  a  3  x  3  second  order  tensor  and  can  be  easily  inverted  as  long  as  the  elements 
are  convex.  The  error  is  defined  as 

error  =  ||/(£")||2  =  )  ■  (12) 

For  the  given  problem,  this  method  typically  converges  to  an  error  of  less  than  1  x 
10“ 6  within  1-3  iterations. 


The  coordinates  of  the  vascular  point  in  the  element  reference  domain  can  then  be 
compared  with  the  bounds  of  the  element.  For  the  hexahedral  element,  the 
coordinates  must  fall  between  -1  and  1 : 

-l<a,77,0<l.  (13) 

And  for  the  tetrahedral  element  the  coordinates  must  fall  between  0  and  1 : 

0<a,77,O<l.  (14) 

If  the  coordinates  fall  within  that  range,  then  the  point  is  contained  in  the  current 
element,  and  the  search  can  tenninate.  If  any  of  the  coordinates  for  the  vascular 
point  are  outside  of  that  range,  the  point  is  not  within  the  current  element,  and  the 
search  algorithm  moves  on  to  the  next  element  and  continues  with  the  coarse 
search. 


2.3  Preprocessing  Algorithm 

In  addition  to  the  goal  of  linking  the  solid  element  head/brain  model  with  the  beam 
element  vascular  model  to  predict  the  strains  and  potential  damage  to  the 
vasculature,  the  vascular  model  can  also  be  used  to  improve  the  accuracy  of  the 
head/brain  model.  The  vascular  structure  within  the  brain  and  along  its  surface 
provides  a  stiffening  effect  and  will  influence  its  defonnation.  In  most  numerical 
models  the  vascular  structure  is  not  included.  Instead  the  brain  matter  and 
vasculature  are  merged  into  a  single  homogenous  material  with  properties  stiffer 
than  white  or  grey  matter  would  have  on  their  own.  Although  this  approximation  is 
commonly  used,  it  does  not  take  into  account  the  structured  nature  of  the 
vasculature,  and  the  effect  is  applied  as  an  isotropic  increase  in  stiffness  instead  of 
the  more  appropriate  anisotropic  increase  in  stiffness.  The  presence  of  the 
vasculature  can  be  included  as  an  additional  anisotropic  tenn  in  the  constitutive 
model  for  the  brain  tissue  based  on  the  direction  of  the  vessels  embedded  within  a 
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given  element.  However,  let  us  begin  by  looking  at  the  material  model  for  the  brain 
prior  to  the  inclusion  of  the  vasculature  tenn. 

In  many  biological  tissues,  fibers  or  bundles  of  cells  are  aligned  in  uniform 
directions.  As  a  result,  isotropic  materials  models  are  insufficient  for  capturing  the 
mechanical  behavior.  For  the  white  matter  within  the  brain,  axon  bundles  form 
complex  fiber  tracts  as  they  connect  and  facilitate  communicate  between  different 
regions  in  the  brain.  These  axonal  fiber  tracts  have  been  reported  to  be 
approximately  3  times  stiffer  than  the  surrounding  matrix  material  and  thus  play  an 
important  role  in  the  mechanical  response  of  the  brain  (Arbogast  and  Margulies 
1999).  To  model  the  white  matter  we  will  use  a  transversely  isotropic  hyperelastic 
material,  where  the  fiber  tract  directions  are  determined  from  diffusion  tensor 
imaging  (DTI)  data  (Kraft  and  Dagro  2011). 


To  describe  the  material  model  we  must  first  define  some  basic  kinematic  concepts. 
The  deformation  gradient  is  defined  as 


F 


dx 
dX  ’ 


(15) 


where  X_  is  the  position  of  a  material  point  in  the  reference  (undefonned) 
configuration  and  x  is  the  position  of  the  same  material  point  in  the  current 
(defonned)  configuration.  The  reference  configuration  refers  to  the  undefonned 
physical  domain,  not  the  element’s  reference  domain.  The  ratio  of  the  deformed 
volume  to  the  undefonned  volume  is  given  by  the  Jacobian,  the  determinant  of  the 
defonnation  gradient 

]  =  det(F (16) 

It  is  often  beneficial  to  perfonn  a  multiplicative  decomposition  of  F  into  volume¬ 
changing  (dilatational)  and  volume-preserving  (distortional)  parts  to  separate  the 
bulk  and  the  shear  response.  To  accomplish  this,  a  deviatoric  deformation  gradient 
in  which  the  volume  change  is  eliminated  is  defined  as 

f=j-l/3p  (17) 


We  can  then  define  a  modified  right  Cauchy-Green  tensor 

C  =  FF.  (18) 


The  modified  principle  invariants  of  the  right  Cauchy-Green  deformation  tensor  are 
defined  as 

=  trC_ , 
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(19) 


(20) 


'2=j[(tr£)2-tr(£2)  , 
and 

J3  =  detC_  =  [detF^j  =  1  .  (21) 

For  the  modified  right  Cauchy-Green  tensor,  because  the  volume  change  has  been 
eliminated,  the  third  principle  invariant  will  always  be  one.  To  capture  the 
anisotropic  nature  of  the  white  matter  we  introduce  a  unit  vector  a0,  assigned  using 
DTI  data  that  describes  the  direction  of  the  fiber  in  the  undefonned  reference 
configuration.  We  can  then  define  2  additional  invariants  based  on  the  fiber 
direction 

I4  —  a0  '  C_a0  ,  (22) 

h  =  9lo  ‘  £2«o  ,  (23) 

where  74  and  75  arise  from  the  anisotropy  and  describe  the  deformation  of  the  fiber 
family.  It  should  be  noted  that 

J4  =  a0-  Ca0  =  J~2/3a  ■  a  =  J~2/3A2  ,  (24) 

where  a  —  F_a0  is  the  direction  of  the  fiber  in  the  current  configuration  and  A  is  the 

stretch  in  the  fiber  bundle.  Thus,  74  will  also  be  useful  in  evaluating  strain  based 
injury  criteria  for  the  axon  fiber  bundles  (Kleiven  2007).  To  take  into  account  the 
presence  of  the  vasculature  we  define  2  more  invariants 

f 6  =  "  —  —0  >  (25) 

1?  —  b0  ■  C_  b0  ,  (26) 

where  b0  is  the  direction  associated  with  the  vasculature  in  the  undefonned 
configuration. 

The  strain  energy,  T,  of  the  transversely  isotropic  hyperelastic  material  can  be 
written  as  a  function  of  the  modified  principle  invariants  along  with  the  Jacobian, 
which  describes  the  change  in  volume.  Assuming  that  the  responses  of  the  fibers 
and  the  matrix  material  are  not  strongly  coupled,  we  can  choose  to  separate  the 
strain  energy  into  a  linear  combination  of  the  isotropic  and  anisotropic  components 

VQl.hjjJs)  =  Visopr  ,h,j)  +  ^aniso  QJs),  (27) 

where  T'iso(71,72,y)  describes  the  response  of  the  isotropic  matrix  and 
TJaniso(74. 75)  describes  the  directional  contribution  of  the  reinforcing  fiber 
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bundles.  We  can  then  select  an  appropriate  isotropic  strain  energy  function  for  the 
matrix  component  such  as  a  Neo-Hookean  or  Mooney-Rivlin  material  model.  For 
the  anisotropic  response  it  is  suggested  to  select  a  Fung  material  model  that  includes 
the  exponential  behavior  characteristic  of  most  soft  tissues 

^anisoft)  =  K  [ exp  (fc2(J4  -  1))  -  J4],  (28) 

where  /q  and  k2  are  material  constants  obtained  from  a  parameter  fit  to 
experimental  data.  This  is  purely  a  simple  example  of  a  Fung  material.  Depending 
on  the  available  data  a  more  complex  constitutive  model  for  the  fibers  may  be 
preferred  (Weiss  et  al.  1996;  Holzapfel  2000;  Ning  et  al.  2006). 


As  with  the  axonal  fibers,  the  direction  of  the  vessels  embedded  within  a  given 
element  can  be  used  to  add  an  anisotropic  component  to  the  constitutive  equation. 
In  the  case  where  only  a  single  vessel  segment  is  contained  in  the  element,  the 
direction  is  simply  determined  from  the  position  of  the  vessel  segment’s  2  end 
points.  However,  there  will  often  be  cases  where  multiple  segments  are  contained 
within  a  single  element.  In  those  cases  a  single  average  direction  must  be 
detennined.  This  can  be  done  by  adding  the  directions  of  the  segments  and  then 
nonnalizing  to  create  a  unit  vector.  However,  the  ordering  of  the  nodes  within  a 
segment  can  cause  the  direction  of  the  /th  segment  to  be  either  dL  or  -  since 


dt 


(29) 


where  x[  and  x_z  are  the  locations  of  the  nodes  on  either  end  of  the /th  segment.  Two 
parallel  segments  of  a  vessel  could  potentially  cancel  each  other  out  depending  on 
the  number  of  the  nodes  (Fig.  5).  To  prohibit  this  from  happening,  an  additional 
tenn  is  included  in  the  vector  sum  to  flip  the  directions  as  needed.  The  average 
vessel  direction  is  calculated  with  the  following  equation: 


H?=l  sid_i 
\'^i=lsid_i 


(30) 


where  sL  is  a  sign  tenn  (sL  =  ±  1)  described  by  the  equation 

Sj  =  1;  i  —  1 

-  1  (31) 

\diilp1sjdj)  I’ 


For  the  initial  segment  s{  —  1,  and  for  all  the  following  segments  st  is  used  to 
detennine  whether  or  not  to  flip  the  /th  segment  based  on  the  sum  of  the  preceding 
segment  directions. 
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Fig.  5  Depending  on  order  of  the  nodes  in  each  vessel  element,  2  parallel  vessel  segments 
could  have  opposite  directions.  In  the  vector  sum  these  will  cancel,  so  the  direction  segment 
must  therefore  be  flipped. 


Using  the  average  vessel  direction,  we  can  calculate  our  6th  invariant  (/6  = 
b0  ■  C_  b0)  and  use  it  to  create  an  additional  term  within  the  anisotropic  portion  on 

the  constitutive  equation.  As  with  the  response  of  the  fiber  tracts,  we  employ  a 
simple  exponential  Fung  type  material  model 


'PanisoO^ 4)  =  K  [exp  (fc2(/4  -  1))  -  /4]  +  k3  [ exp  (fc4(/6  -  l))  -  /6].  (32) 


The  calculation  for  the  average  direction  works  well  when  summing  multiple 
segments  that  are  close  to  parallel.  However,  when  segments  are  nearly 
perpendicular  an  average  direction  does  not  have  much  physical  significance.  Take 
the  example  of  an  element  containing  3  perpendicular  fiber  segments  of  equal 
length.  For  the  sake  of  simplicity,  let  those  sections  be  aligned  with  the  x,  y,  and  z 
axes.  The  average  direction  would  be  the  sum  of  these  vectors: 


§_x  @y  T  (Z7 


T 

0 

+ 

O' 

1 

+ 

O' 

0 

= 

T 

1 

.0. 

.0. 

.1. 

.1. 

(33) 


And  that  sum  is  then  normalized  to  calculate  the  directional  vector  associated  with 
that  element: 


(34) 


However,  because  these  vectors  making  up  the  sum  are  perpendicular,  other  vectors 
would  be  equally  valid  as  the  segment  vectors  could  be  defined  as  ±_e,  and  thus 


‘o  =  ±Vf 


-r 

- 1 ' 

■  l ' 

l 

-l 

ii 

i+ 

U- 

l 

.  l . 

.  l . 

.-l. 

(35) 


Approved  for  public  release;  distribution  is  unlimited. 

12 


are  also  possible  solutions.  In  truth,  there  should  be  no  preferred  direction  for  this 
element.  So  we  introduce  a  scale  factor  to  allow  the  user  to  take  this  into  account 
when  creating  the  material  model.  This  scaling  factor  is  defined  as  the  sum  of  the 
absolute  values  of  the  dot  product  of  the  segments,  dh  with  the  previously 
calculated  directional  vector,  h0,  divided  by  the  sum  of  the  lengths  of  the  segment 
vectors 


f  _  2f=ilko -di 

1  ~  ZIUI&II 


(36) 


This  allows  the  length  of  each  directional  vector  to  act  as  a  weight  in  the  sum.  The 
maximum  value  for  the  scaling  factor,  when  all  of  the  segments  are  aligned,  is 
therefore  one.  For  the  special  case  where  there  are  3  segment  vectors  of  equal  length 

but  perpendicular  to  one  another,  the  scale  factor  would  be  or  0.57735.  As  a 

1 

result  we  expect  our  scaling  factor  to  have  a  range  of  <  /  <  1.  We  can  adjust 
that  range  by  defining  a  new  scaling  factor 


c  3+77  fgjLiferTil  _  j_\ 

'  2  UlUNI  V3 ) 


(37) 


such  that  0  <  /  <  1.  This  scaling  factor  is  a  measure  of  the  correlation  between 
fiber  directions  and  can  be  considered  a  measure  of  the  anisotropy  introduced  by 
the  presence  of  the  vasculature.  Highlight-aligned  vasculature  will  leave  to  a 
scaling  factor  of  close  to  1  and  have  a  highly  anisotropic  effect  on  the  overall 
response.  Conversely,  evenly  distributed  perpendicular  vessels  will  lead  to  a  scaling 
factor  of  close  to  zero  and  have  a  more  isotropic  effect  on  the  overall  response.  The 
preprocessing  algorithm  reports  the  scaling  factor  for  each  element  along  with  the 
directional  vector,  b_Q.  It  is  suggested  that  the  updated  scaling  factor,  /,  be  used  to 
weight  the  anisotropic  portion  of  the  constitutive  equation  that  corresponds  to  the 
vascular  structure 


^amso(U.U)  =  K  [exp  (k2(! 4  -  1))  -  z4]  +  fh  [exp  ( k4(l6  -  l))  -  /6]  . 

(38) 


To  account  for  the  isotropic  stiffening  effect  of  perpendicularly  aligned  vessels 
(cases  where  /  is  close  to  zero),  the  isotropic  component  of  constitutive  equation 
can  also  be  updated 

'VQlJlJ.UJs)  =  [1  +  K*(l  -/)]^iso(7i,72j)  T^niso^/s)  •  (39) 

Here  the  material  constant,  k*,  represents  the  traditional  artificial  isotropic 
stiffening  of  the  brain  tissue  due  to  the  presence  of  the  vessels.  As  the  vessels  are 
more  unifonnly  aligned,  that  term  goes  to  zero  and  the  anisotropic  term  takes  over. 
If  the  vessels  are  randomly  oriented  and  /  approaches  zero,  the  additional 
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anisotropic  term  vanishes  and  isotropic  stiffness  is  increased.  This  code  merely 
provides  a  calculation  for  the  average  vessel  direction  and  the  scale  factor.  The 
decision  on  whether  or  not  to  use  this  data  to  alter  the  material  model  and  how  it  is 
used  within  the  constitutive  equation,  is  left  to  the  end  user. 


2.4  Postprocessing  Algorithm 


To  calculate  the  potential  damage  in  the  vascular  structure  we  must  be  able  to  map 
to  the  defonnation  from  the  results  of  the  brain  simulation  onto  the  vascular  mesh. 
The  algorithm  described  in  Section  2.2  has  determined  the  solid  element  containing 
each  vessel  segment.  It  has  also  detennined  the  corresponding  coordinates  in  the 
element’s  reference  domain.  Combining  this  information  with  the  output  data  from 
the  simulation  of  the  brain  simulation  will  allow  us  to  determine  the  strains  within 
the  vascular  network.  Although  the  initial  locations  of  the  vascular  nodes  are  not 
related  to  the  nodes  in  the  brain  mesh,  the  vascular  nodes  are  treated  as  material 
points  within  the  brain  mesh.  Thus,  as  the  brain  mesh  deforms,  so  too  do  the 
vascular  nodes.  It  is  that  deformation  which  the  postprocessing  algorithm  will 
detennine. 


The  finite  element  simulation  produces  an  approximate  solution,  uh,  for  the 
displacement,  which  is  a  function  of  the  degrees  of  freedom,  ul,  and  the  shape 
functions,  Nt.  Both  of  these  correspond  to  specific  nodes  within  the  descritized 
finite  element  mesh. 


NEN 

u(X)  ~  uh{X_ )  =  ^  w'tyQO 

t  =  l 


(40) 


While  the  displacement  is  a  function  of  the  position  within  the  mesh,  the  ul,s  are 
constants.  For  a  given  element,  only  the  shape  functions  corresponding  to  nodes 
within  that  element  have  nonzero  values.  So,  to  detennine  the  deformation  we  only 
need  to  sum  over  the  ul  and  values  conesponding  to  nodes  within  that  element. 
To  determine  the  gradient  of  the  displacement  we  take  advantage  of  the  fact  that 
the  ul,s  are  constants  and  move  the  -j-  within  the  sum.  Thus  the  gradient  of  the 

cix. 

displacement  is  just  the  sum  of  those  constants  times  the  gradient  of  the 
corresponding  shape  function  at  the  current  coordinate 

NEN  ,  . 

du  XT'  ^idNi(X_) 

dX~  2 j-  dX  .... 

-  i=i  -  .  1411 


However,  in  practice  the  shape  functions  are  written  in  terms  of  the  element 
reference  domain’s  coordinates,  so  and  not  the  physical  domain’s 

coordinates,  IVjQf).  The  gradient  of  the  shape  functions  with  respect  to  the  physical 
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coordinates  can  be  rewritten  as  the  gradient  with  respect  to  the  reference  domain’s 
coordinates  multiplied  by  the  derivative  of  the  reference  coordinate  with  respect  to 
the  physical  coordinate 


dNt(x)  _  dW‘(i)  di 
dX  ~  df  dX  ' 


(42) 


For  general  trilinear  hexahedral  elements  it  is  not  possible  to  write  a  simple 
expression  for  in  tenns  of  X;  however,  there  is  a  simple  expression  for  X_  as  a 

function  of 


NEN 

x  =  £  m  @ 

i= 1 


(43) 


where  are  the  same  shape  functions  as  before  and  X_l  are  the  coordinates  of  the 
element’s  nodes  in  the  physical  domain.  This  can  then  be  easily  differentiated  to 
give  us 


dX 

dj 


dN i  (f) 

d<f 


(44) 


The  resulting  3><3  tensor  can  then  be  inverted  to  provide  the  gradient  of  with 
respect  to  X_ 


d£ 

dX 


(45) 


With  the  gradient  of  the  descritized  displacement  determined,  we  can  calculate  the 
strain  at  the  desired  location  in  the  brain  mesh.  This  is  the  strain  that  the 
corresponding  point  in  the  vascular  mesh  will  experience 


1  f  du  du 1 


F  =  - ( —  +  —  +^LT 

2  \dX  dX  dX  dxj 


(46) 


This  calculation  is  for  the  Green-Lagrange  strain  tensor.  While  it  is  beneficial  to 
have  the  entire  strain  tensor,  we  are  especially  interested  in  the  stretch  along  the 
fiber  axis.  As  mentioned  in  Section  2.3,  the  stretch  along  the  fiber  axis  is  related  to 
our  6th  principal  invariant 

I6  =  b0-gb0  =  b-b  =  A2,  (47) 

where  b0  is  the  unit  vector  describing  the  direction  of  the  fiber  and  C_  is  the  right 

Cauchy-Green  tensor.  The  right  Cauchy-Green  tensor  is  related  to  the  strain 
through  the  following  equation: 
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(48) 


£  =  2  E  +  I_. 

The  postprocessing  algorithm  creates  a  data  file  containing  the  element  IDs  for  the 
vessel  segments,  the  6  unique  components  of  the  Green-Lagrange  strain  (Exx,  Eyy, 
Ezz,  Exy,  Eyz,  Ezx),  and  the  stretch  (A)  along  the  segment  axis.  These  data  are 
repeated  for  each  time  step  reported  by  the  output  of  the  solid  mesh  simulation. 

3.  Example 

To  demonstrate  the  capabilities  of  this  code,  a  sample  mesh  was  generated  using 
trilinear  hexahedral  elements.  This  mesh  is  intended  to  represent  a  small  section  of 
the  larger  brain  mesh.  Figure  6  shows  2  views  of  this  mesh.  The  mapping  from  the 
physical  domain  to  the  element’s  reference  domain  is  greatly  simplified  in  the  event 
of  cubic  elements  or  rectangular  cuboids.  Therefore,  the  shape  of  the  mesh  was 
chosen  so  as  to  force  the  individual  elements  to  avoid  these  simpler  shapes  and 
better  demonstrate  the  abilities  of  the  code. 
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Fig.  6  Solid  element  mesh  representing  a  portion  of  the  larger  mesh 

The  vascular  network  mesh  can  be  seen  in  Fig.  7.  This  mesh  occupies  the  same 
physical  space  as  the  solid  element  brain  mesh  but  does  not  share  any  of  its 
geometry.  The  nodes  and  edges  of  the  2  meshes  do  not  necessarily  coincide.  As  the 
vessels  move  away  from  their  starting  point,  they  branch  and  alter  direction, 
creating  irregular  paths  through  space. 
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Fig.  7  Vascular  structure  made  up  of  1-D  elements 

The  mapping  algorithm  identifies  the  elements  containing  each  of  the  vessel 
segments  and  determines  the  location  of  the  vessel  in  terms  of  the  element’s 
reference  coordinates.  Figure  8  shows  the  vascular  network  overlaid  on  the  solid 
element  mesh.  Using  the  results  of  the  mapping  algorithm,  only  the  elements 
containing  a  vessel  segment  are  displayed.  For  this  example  the  vessels  continue 
on  their  paths  until  they  have  left  the  solid  mesh.  Endpoints  of  the  vessels  in  the 
image  on  the  right  of  Fig.  8  can  be  seen  outside  of  any  element.  This  is  not  an  error 
in  the  code,  this  is  because  they  are  outside  the  boundaries  of  the  solid  mesh  and 
thus  are  not  contained  within  any  element.  These  free  endpoints  will  not  be  a  factor 
in  any  of  the  calculations  in  either  the  pre-  or  postprocessing  algorithm. 
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Fig.  8  Vascular  structure  and  solid  elements  containing  vessels 


3.1  Preprocessing  Results 

Using  the  preprocessing  algorithm,  every  element  containing  a  vessel  segment  was 
assigned  a  directional  vector,  b0,  and  a  scale  factor  as  a  measure  of  the  correlation 
between  the  directions  of  the  vessel  segments.  Table  2  shows  a  sample  of  the  data 
reported  by  the  preprocessing  algorithm.  Because  of  the  irregular  paths  and 
branching  vessels,  the  scale  factor  shows  a  value  less  than  one  when  multiple 
segments  were  present  within  a  single  element  due  to  variations  in  the  directions  of 
the  individual  segments. 
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Table  2  Preprocessing  data  for  a  sample  of  the  solid  elements  containing  vascular  segments 


Element 

X 

Y 

Z 

Scale,  / 

84 

0.369179 

0.897098 

0.771899 

1.0 

156 

0.585461 

0.375538 

0.718475 

0.549186 

157 

0.455495 

0.470486 

0.755756 

0.791908 

164 

0.213751 

0.930328 

0.297994 

0.819901 

172 

0.737664 

0.670064 

-0.082860 

1.0 

173 

0.933215 

0.340099 

-0.115940 

0.761922 

174 

0.556493 

0.665586 

0.497304 

0.812248 

182 

0.303811 

0.851047 

0.428273 

0.389745 

190 

0.554413 

0.826289 

-0.099359 

1.0 

222 

0.793040 

0.311000 

0.523800 

1.0 

230 

0.193325 

0.740923 

0.643163 

0.573547 

231 

0.581279 

0.548160 

0.601362 

0.753829 

3.2  Postprocessing  Results 

To  examine  the  mapping  of  strains  from  the  solid  mesh  to  the  vessel  network,  the 
solid  mesh  was  subjected  to  a  time  varying  strain  that  also  varied  nonlinearly  with 
position.  The  displacement  values  are  defined  at  the  nodes  of  the  solid  mesh  only, 
not  the  nodes  of  the  vascular  mesh.  The  displacement  field  defined  over  the  solid 
mesh  was  then  mapped  onto  the  vascular  mesh,  and  the  strains  in  the  vessel 
segments  were  calculated.  Figure  9  shows  a  plot  of  the  vascular  network  with  the 
color  of  the  vessel  segments  indicating  the  amount  of  stretch  in  each.  As  expected, 
the  vessel  segments  in  the  regions  where  solid  elements  have  the  largest  strain  also 
exhibited  the  largest  strain  values.  However,  fiber  direction  plays  a  role  as  well 
since  the  stretch  along  the  fiber  axis  is  a  function  of  the  fiber’s  direction 

A  =  Jb0-£b0  .  (49) 

Thus,  there  is  additional  variation  in  the  stretch  and  axial  strain  within  the  fibers 
due  to  variations  in  their  orientations.  Merely  calculating  the  strain  within  the  solid 
mesh  does  not  provide  enough  information  to  make  an  accurate  prediction  of  the 
potential  for  damage  in  the  vasculature;  the  direction  of  the  vessel  segments  must 
be  taken  into  account. 
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0.8  _ 


Fig.  9  Axial  stretch  in  the  vessels,  ranging  from  1  (undeformed)  to  1.5  (50%  elongation) 

4.  Conclusion 


TBI  is  a  serious  concern  in  both  the  military  and  civilian  population.  Even  when 
nonfatal,  injury  to  the  brain  has  the  potential  for  lifelong  repercussions.  For  the 
Soldier,  blast-induced  TBI  is  of  particular  concern.  While  significant  research  has 
been  done  on  models  to  predict  TBI,  many  of  them  have  focused  on  axonal  injury 
and  often  neglect  the  possible  presence  of  damage  to  the  cerebral  vasculature. 
During  Operation  Iraqi  Freedom,  a  study  of  Soldiers  at  Walter  Reed  Army  Medical 
Center  suffering  TBI  found  that  47%  had  vasospasm,  35%  had  a  pseudoaneurysm, 
12%  had  a  subarachnoid  hemorrhage,  3.5%  had  an  epidural  hematoma,  16%  had  a 
subdural  hematoma,  11%  had  an  intraventricular  hemorrhage,  and  14%  had  mixed 
hemorrhages  (Armonda  et  al.  2006).  Injuries  such  as  vasospasm  and 
pseudoaneurysm  can  lead  to  further  damage  to  the  brain  over  time.  Intracranial 
hemorrhages  and  hematomas  can  be  life  threatening  and  often  require  immediate 


Approved  for  public  release;  distribution  is  unlimited. 

21 


intervention.  While  damage  to  the  vasculature  is  a  common  feature  of  TBI,  many 
numerical  models  for  the  brain  do  not  include  the  vascular  structures  within  the 
brain  and  thus  are  incapable  of  predicting  damage  to  the  cerebral  vasculature.  This 
work  developed  a  code  capable  of  postprocessing  an  existing  head  and  brain  model 
to  detennine  the  deformation  in  a  corresponding  cerebral  vasculature  mesh.  With 
this  infonnation,  the  user  can  then  input  the  stretches  and  strains  of  the  vasculature 
into  a  damage  model  of  their  choice  to  predict  potential  vascular  injury.  The  code 
is  also  capable  of  using  the  geometry  of  the  cerebral  vasculature  to  provide 
information  for  an  updated  anisotropic  material  model,  allowing  the  directional 
nature  of  the  vessels  to  inform  the  material  response  of  the  surrounding  brain  tissue. 

This  approach  is  particularly  desirable  because  it  builds  upon  an  existing  head  and 
brain  model  currently  in  use  by  the  Soldier  Protection  Sciences  Branch  at  the  US 
Anny  Research  Laboratory.  Whereas  developing  an  entirely  new  head  model  and 
mesh  would  be  extremely  time  consuming,  this  code  can  be  immediately  put  into 
use  with  the  existing  simulations.  Another  important  feature  of  this  model  is  that 
the  topography  of  the  brain  and  vascular  meshes  can  be  independent  of  one  another. 
Other  approaches  to  linking  a  solid  mesh  with  1 -dimensional  (1-D)  network  often 
require  that  the  locations  of  the  nodes  and  edges  in  the  2  meshes  coincide.  For  the 
complex  geometry  of  the  brain  and  cerebral  vasculature,  this  requirement  would 
impose  severe  restrictions  on  the  brain  mesh,  potentially  increasing  simulation  time 
or  even  making  effective  meshing  impossible.  This  approach  places  no  such 
restrictions  on  either  mesh  and  requires  no  changes  to  the  existing  head  and  brain 
model.  Another  important  feature  of  this  code  is  that  it  can  be  used  to  help  improve 
the  material  model  for  the  brain.  The  presence  of  the  vasculature  within  the  brain 
produces  a  reinforcing  effect  and  leads  to  an  increase  in  the  effective  stiffness  of 
the  combined  brain-vascular  material.  While  many  simulations  account  for  this 
with  isotropic  increases  to  the  stiffness  of  the  brain’s  material  model,  this  code  can 
be  used  to  assign  directional  values  to  that  increased  stiffness  and  lead  to  a  more 
accurate  anisotropic  material  model.  Finally,  this  code  gives  the  user  a  great  deal 
of  flexibility  in  how  they  determine  the  damage  in  the  vasculature.  Since  the  raw 
strain  and  axial  stretch  data  for  the  vessel  segments  are  reported  at  each  output  step, 
the  user  can  use  whatever  damage  model  they  deem  most  appropriate.  Although 
this  approach  has  many  strengths,  there  are  some  limitations.  The  model  does  not 
currently  account  for  fluid  flow  within  the  arteries  and  the  corresponding  internal 
pressure  of  the  vessels.  Also,  since  it  uses  a  1-D  approximation  of  the  vessels,  any 
damage  model  must  be  simplified.  However,  both  of  these  limitations  can  be 
addressed  through  future  work. 

There  are  a  number  of  new  features  that  would  be  desirable  if  work  on  the  project 
continues.  In  regards  to  the  preprocessing  algorithm,  inclusion  of  data  on  the 
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thickness  of  the  vessels  could  be  used  to  provide  more  accurate  anisotropic  material 
constants.  Larger  vessels  would  have  a  more  pronounced  effect  on  the  stiffening  of 
the  surrounding  brain  tissue;  therefore  the  thickness  of  the  segments  could  be  used 
to  provide  an  additional  scaling  tenn  for  the  constitutive  equation.  Furthermore,  a 
secondary  simulation  could  be  created  to  model  the  fluid  flow  within  the  arteries. 
Using  a  pressure  gradient  to  drive  the  blood  flow,  and  the  external  pressure  induced 
by  a  blast  wave  through  the  surrounding  brain  elements,  an  approximation  for  the 
time  varying  internal  pressure  of  the  vasculature  could  theoretically  be  created. 
When  this  internal  pressure  calculation  is  used  along  with  a  mapping  of  the  external 
pressure  from  the  brain,  a  measure  of  the  compressive  stress  through  the  thickness 
of  the  vascular  walls  could  be  attained.  This  compressive  wall  stress  is  likely  an 
important  contributor  to  and  predictor  of  vasospasm  in  the  arteries.  Finally,  future 
work  could  directly  model  the  largest  cerebral  vascular  structures  in  the  head  and 
brain  model,  relying  on  this  approach  for  only  the  smaller  scale  vasculature. 
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To  read  the  mesh  files  used  in  SIERRA,  this  program  requires  the  EXODUS  II  C 
and  C++  libraries.  The  files  necessary  to  generate  the  EXODUS  II  library  can  be 
found  at  sourceforge.net.  These  must  be  compiled  and  the  location  of  the  EXODUS 
library  must  be  linked  when  the  cerebral  vascular  injury  code  is  compiled.  It  is  also 
possible  to  read  in  the  mesh  data  using  LS-DYNA  keyword  format;  however,  the 
EXODUS  II  library  is  still  necessary  to  compile  the  code.  If  the  user  wishes  to  avoid 
using  the  Exodus  libraries,  it  is  recommended  to  create  a  separate  copy  of  this  code 
with  the  “include<exodusII.h>”  line  removed  or  commented  out  and  all  of  the 
EXODUS  specific  functions  in  the  body  of  the  code. 

A-l  General  Input  File  Commands 

$  at  the  beginning  of  a  line  denotes  a  comment.  Any  infonnation  following  the  $ 
will  be  ignored 

*PREPROCESS  ON/OFF 

Detennines  whether  the  preprocessing  algorithm  is  run  to  assign  average  directions 
to  the  brain  elements.  Acceptable  inputs  are  on,  On,  ON,  off,  Off,  or  OFF.  If  the 
command  is  omitted,  the  simulation  treats  it  as  off. 

*POSTPROCESS  ON/OFF  <fonnat_type>  filename,  ex  tension 
Detennines  whether  the  postprocessing  algorithm  is  run  to  calculate  the  strain  in 
the  vascular  beam  elements.  Acceptable  inputs  are  on,  On,  ON,  off,  Off,  or  OFF.  If 
the  command  is  omitted,  the  simulation  treats  it  as  off.  After  the  command  is  the 
format  type  for  the  data  and  the  name  of  the  file  containing  the  finite  element 
deformation  data  for  the  solid  mesh.  The  fonnat  type  is  an  integer  value  and  must 
be  included.  The  acceptable  values  are  as  follows: 

1  -  For  output  in  EXODUS  II  format 

2  -  For  output  in  text  format  (for  more  information  see  Appendix  C) 

The  file  name  can  include  the  path  to  the  file  if  it  is  not  in  the  current  directory. 

*NEWTON 
$  max_it  tol 
20  le-6 

Sets  parameters  for  Newton  search  algorithm. 

Max_it  -  Maximum  number  of  Newton  iterations  (integer,  default  =  20) 
tol  -  Tolerance  for  Newton  iteration  (double,  default  le-6) 

*INCLUDE  file  name. extension 
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Allows  for  the  input  file  to  include  other  files.  Primarily  used  so  that  the  mesh  can 
be  stored  in  a  separate  file.  Multiple  files  can  be  included  and  those  files  may  also 
contain  the  include  command.  The  name  of  the  file  and  its  extension  must  follow 
the  *  INCLUDE  command 

*END 

Denotes  the  end  of  an  input  file.  Can  be  either  the  end  of  the  main  input  file  or  of 
an  include  file. 

A-2  Mesh  Input  File  Commands 

It  is  possible  to  read  in  mesh  data  either  through  importing  a  file  in  EXODUS  II 
format  or  using  LS-DYNA  keyword  input  fonnat. 

* EXODU S  INPUT  S OLID  file  name .  extension 

Read  in  the  mesh  file  containing  the  solid  elements  for  the  brain/head  model.  File 
must  be  in  EXODUS  II  format. 

*  EXODUS_INPUT_BEAM  filename,  ex  tension 

Read  in  the  mesh  file  containing  the  beam  elements  for  the  vasculature  model.  File 
must  be  in  EXODUS  II  format. 

Keyword  input  files  may  be  used  in  place  of  the  EXODUS  II  fonnat  files.  The 
following  are  the  commands  to  input  the  nodes  and  elements  for  either  mesh.  The 
simulation  will  function  with  either  fonnat,  you  do  not  need  to  use  both. 

*NODE 


$  NID 

X 

Y 

1 

0.0 

0.0 

0.0 

2 

1.0 

0.0 

0.0 

3 

0.0 

1.0 

0.0 

4 

1.0 

1.0 

0.0 

Etc. 


NID  -  Node  ID  Number  (integer  >  0):  Unique  identifier  for  each  node.  Numbers 
must  be  greater  than  zero,  however  numbering  does  not  need  to  begin  at  one  or  be 
sequential. 

X,  Y,  Z  -  Reference  position  for  the  node  (double) 

*  ELEMENT  S  OLID 

$  EID  type  NIDI  NID2  NID3  NID4  NID5  NID6  NID7  NID8 

1  1  1  2  12  11  41  42  48  47 
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2  1  2  3  13  12  42  43  49  48 

Etc. 

EID  -  Element  ID  Number  (integer  >  0):  Unique  identifier  for  each  element. 
Numbers  must  be  greater  than  zero,  however  numbering  does  not  need  to  begin  at 
one  or  be  sequential 

Type  -  Element  type  (integer):  Used  in  LS-DYNA  but  not  only  a  placeholder  here. 
Some  value  must  be  entered  but  actual  value  is  irrelevant 

NIDI  -  NID8  -  Node  ID  Numbers  corresponding  to  the  vertices  of  the  element.  A 
hexahedral  element  will  have  8  NIDs  defined.  A  tetrahedral  element  will  have  4 
NIDs  defined.  Based  on  the  number  of  NIDs  entered,  the  simulation  will 
automatically  detennine  if  a  given  solid  element  is  a  hexahedral  or  tetrahedral 
element.  All  NID  numbers  must  correspond  to  nodes  defined  using  *NODE 

*ELEMENT_BEAM 

$  Beam  Elements  that  make  up  the  vascular  network 
$  EID  type  NIDI  NID2 
1112 
2  13  4 

Etc. 

EID  -  Element  ID  Number  (integer  >  0):  Unique  identifier  for  each  element. 
Numbers  must  be  greater  than  zero,  however  numbering  does  not  need  to  begin  at 
one  or  be  sequential.  Ideally  beam  elements  would  not  share  EIDs  with  solid 
elements;  however,  if  they  do  it  should  not  cause  a  problem  with  the  simulation. 
Type  -  Element  type  (integer):  Used  in  LS-DYNA  but  not  only  a  placeholder  here. 
Some  value  must  be  entered  but  actual  value  is  irrelevant 

NIDI  -  NID2  -  Node  ID  Numbers  corresponding  to  the  vertices  of  the  element.  All 
NID  numbers  must  correspond  to  nodes  defined  using  *NODE 
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Appendix  B.  Sample  Input  File 
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B-l  Primary  Input  File 


*PRE PROCESS  ON 
$ 

*POSTPROCESS  ON  2  post_data . output 
$ 

*NEWTON 
$  maxjt  tol 
20  le-6 

*  INCLUDE  Mesh. input 
*END 


B-2  Include  File  with  Mesh 


Input  file  using  LS-DYNA  keyword  fonnat  for  the  mesh.  EXODUS  II  input  file 
can  be  used  in  lieu  of  keyword  input  (see  Appendix  A  Section  A-2). 

<Mesh.input> 

*NODE 

5  0.250000  0.250000  0.000000 

6  0.500000  0.125000  0.000000 

7  0.750000  0.000000  0.000000 

8  0.250000  0.375000  0.000000 

9  0.562500  0.437500  0.000000 

10  0.875000  0.500000  0.000000 

11  0.250000  0.500000  0.000000 

12  0.625000  0.750000  0.000000 

13  1.000000  1.000000  0.000000 

14  0.250000  0.250000  0.500000 

15  0.500000  0.125000  0.500000 

16  0.750000  0.000000  0.500000 

17  0.250000  0.375000  0.500000 

18  0.562500  0.437500  0.500000 

19  0.875000  0.500000  0.500000 

20  0.250000  0.500000  0.500000 

21  0.625000  0.750000  0.500000 

22  1.000000  1.000000  0.500000 

23  0.250000  0.250000  1.000000 

24  0.500000  0.125000  1.000000 

25  0.750000  0.000000  1.000000 

26  0.250000  0.375000  1.000000 

27  0.562500  0.437500  1.000000 

28  0.875000  0.500000  1.000000 

29  0.250000  0.500000  1.000000 
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30  0.625000  0.750000  1.000000 

31  1.000000  1.000000  1.000000 

32  0.273756  0.438769  0.210000 

33  0.391502  0.420442  0.362500 

34  0.509247  0.402116  0.515000 

35  0.626992  0.383789  0.667500 

36  0.744737  0.365463  0.820000 

37  0.273756  0.438769  0.210000 

38  0.321508  0.375461  0.387500 

39  0.369259  0.312153  0.565000 

40  0.417011  0.248845  0.742500 

41  0.464762  0.185537  0.920000 
*ELEMENT_SOLID 

3  1  5  6  9  8  14  15  18  17 

4  1  6  7  10  9  15  16  19  18 

5  1  8  9  12  11  17  18  21  20 

6  1  9  10  13  12  18  19  22  21 

7  1  14  15  18  17  23  24  27  26 

8  1  15  16  19  18  24  25  28  27 

9  1  17  18  21  20  26  27  30  29 

10  1  18  19  22  21  27  28  31  30 
*ELEMENT_BEAM 

11  1  32  33 

12  1  33  34 

13  1  34  35 

14  1  35  36 

15  1  37  38 

16  1  38  39 

17  1  39  40 

18  1  40  41 
*END 

B-3  Post-Process  Data  File 


<p°st_data.output> 


3 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 


0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 


0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 
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0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.100000 

0.012500 

0.006250 

0.000000 

0.018750 

0.021875 

0.025000 

0.025000 

0.037500 

0.050000 

0.012500 

0.006250 

0.000000 

0.018750 

0.021875 

0.025000 

0.025000 


0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.025000 

0.012500 

0.000000 

0.037500 

0.043750 

0.050000 

0.050000 

0.075000 

0.100000 

0.025000 

0.012500 

0.000000 

0.037500 

0.043750 

0.050000 

0.050000 


0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.050000 

0.050000 

0.050000 

0.050000 

0.050000 

0.050000 

0.050000 
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0.037500 

0.050000 

0.012500 

0.006250 

0.000000 

0.018750 

0.021875 

0.025000 

0.025000 

0.037500 

0.050000 

0.021938 

0.021022 

0.020106 

0.019189 

0.018273 

0.021938 

0.018773 

0.015608 

0.012442 

0.009277 

0.200000 

0.025000 

0.012500 

0.000000 

0.037500 

0.043750 

0.050000 

0.050000 

0.075000 

0.100000 

0.025000 

0.012500 

0.000000 

0.037500 

0.043750 

0.050000 

0.050000 

0.075000 

0.100000 

0.025000 

0.012500 

0.000000 

0.037500 

0.043750 

0.050000 


0.075000 

0.100000 

0.025000 

0.012500 

0.000000 

0.037500 

0.043750 

0.050000 

0.050000 

0.075000 

0.100000 

0.043877 

0.042044 

0.040212 

0.038379 

0.036546 

0.043877 

0.037546 

0.031215 

0.024885 

0.018554 

0.050000 

0.025000 

0.000000 

0.075000 

0.087500 

0.100000 

0.100000 

0.150000 

0.200000 

0.050000 

0.025000 

0.000000 

0.075000 

0.087500 

0.100000 

0.100000 

0.150000 

0.200000 

0.050000 

0.025000 

0.000000 

0.075000 

0.087500 

0.100000 


0.050000 

0.050000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.021000 

0.036250 

0.051500 

0.066750 

0.082000 

0.021000 

0.038750 

0.056500 

0.074250 

0.092000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.100000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 

0.200000 
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0.050000 

0.075000 

0.100000 

0.043877 

0.042044 

0.040212 

0.038379 

0.036546 

0.043877 

0.037546 

0.031215 

0.024885 

0.018554 


0.100000 

0.150000 

0.200000 

0.087754 

0.084088 

0.080423 

0.076758 

0.073093 

0.087754 

0.075092 

0.062431 

0.049769 

0.037107 


0.200000 

0.200000 

0.200000 

0.042000 

0.072500 

0.103000 

0.133500 

0.164000 

0.042000 

0.077500 

0.113000 

0.148500 

0.184000 
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Appendix  C.  Alternative  Format  for  Postprocessing  Data 
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Instead  of  using  an  EXODUS  II  file  for  the  deformation  data,  a  text  file  may  be 
used.  That  file  must  match  the  following  fonnat  for  the  algorithm  to  function.  The 
file  will  contain  only  numbers,  no  keywords  are  used. 

The  first  line  will  be  a  single  integer  entry  indicating  the  number  of  output  time 
steps.  After  that  the  time  of  the  first  output  step  is  reported.  Then  the  following  lines 
include  the  displacement  values  of  the  nodes  of  the  solid  mesh  in  order.  All  nodes 
must  be  reported  at  each  output  time  step. 

<line  1  -  Number  of  output  time  steps> 

<line  2  -  Time  for  output  step  1> 

<line  3  through  number  of  nodes  +  3  -  X,  Y,  and  Z  displacement  values  for  each 
node> 

<line  number  of  nodes  +  4  -  Time  for  output  step  2> 

<line  number  of  nodes  +  5  through  2*  (number  of  nodes)  +  5  -  X,  Y,  and  Z 
displacement  values  for  each  node> 

This  pattern  continues  until  all  of  the  output  time  steps  have  been  reported 
For  an  example  of  this  fonnat,  see  Appendix  B-3. 
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Appendix  D.  EXODUS  II  Commands 


Approved  for  public  release;  distribution  is  unlimited. 


39 


The  following  header  line  is  required  at  the  top  of  the  code  to  use  any  of  the 
following  EXODUS  II  commands: 

#include  “exodusll.h” 

D-l  Reading  in  Mesh  Data 

Open  EXODUS  II  file 

int  ex_open(path,  mode,  comp_ws,  io_ws,  version); 

path  -  char,  file  name  and  path  to  the  EXODUS  II  file 

mode  -  int,  EX  READ  or  EX  WRITE  (predefined  constants) 

c°mp_ws  -  int,  word  size  in  bytes  (0,  4,  or  8) 

io_ws  -  int,  word  size  in  bytes  (0,  4,  or  8) 

version  -  float*,  returned  EXODUS  II  database  version  number 

Close  EXODUS  II  file 

int  ex_close  (exoid) ; 

exoid  -  int,  returned  value  from  ex  open 

Read  Initialization  Parameters 

int  ex_get_init( exoid,  title,  numdim,  numnodes,  numelem,  numelemblk, 
num_node_sets,  num_side_sets); 

exoid  -  int,  returned  value  from  ex  open 

title  -  char*,  returned  database  title 

num  dim  -  int*,  returned  dimensionality  of  the  database 

num  nodes  -  int*,  returned  number  of  nodes 

num  elem  -  int*,  returned  number  of  elements 

num  elem  blk  -  int*,  returned  number  of  element  blocks 

numnodesets  -  int*,  returned  number  of  node  sets 

num  side  sets  -  int*,  returned  number  of  side  sets 

Read  Nodal  Coordinates 

int  ex_get_coord( exoid,  xcoor,  y  coor,  zcoor); 

exoid  -  int,  returned  value  from  ex  open 
x  coor  -  float*,  returned  x  coordinates  of  the  nodes 
y  coor  -  float*,  returned  y  coordinates  of  the  nodes 
z  coor  -  float*,  returned  z  coordinates  of  the  nodes 
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Read  Element  Block  Parameters 

int  ex_get_elem_block(exoid,  elem_blk_id,  elemtype,  numelemthisblk, 
num _nodes_per_elem,  numattr); 

exoid  -  int,  returned  value  from  ex  open 

elemblkid  -  int,  ID  of  current  element  block 

elem  type  -  char*,  returned  type  of  elements  in  the  element  block 

num  elem  this  blk  -  int*,  returned  number  of  elements  in  the  element  block 

num_nodes_per_elem  -  int*,  returned  number  of  nodes  per  element  in  the  element 

block 

num  attr  -  int*,  returned  number  of  attributes  per  element  in  the  element  block 

Read  Element  Block  Connectivity 

int  ex_get_elem_conn( exoid,  elem  blk  id,  connect); 

exoid  -  int,  returned  value  from  ex  open 
elem  blk  id  -  int,  element  block  ID 

connect[num_elem_this_blk,  num_nodes_per_elem]  -  int,  list  of  nodes  making  of 
the  current  element 

D-2  Reading  in  Result  Data 

Read  Results  Variables  Parameters 

int  ex_get_var_param(exoid,  var_type,  numjvars); 

exoid  -  int,  returned  value  from  ex  open 

var_type  -  char*,  type  of  variable  described 

“g”  or  “G”  -  global  variables 

“n”  or  “N”  -  nodal  variables 

“e”  or  “E”  -  element  variables 

“m”  or  “M”  -  nodeset  variables 

“s”  or  “S”  -  sideset  variables 

numvars  -  int*,  returned  number  of  variables  stored  for  specified  variable  type 

Read  Variable  Names 

int  ex  get  var  names  (exoid,  var  type,  num  vars,  var^names [ ] ) ; 


exoid  -  int,  returned  value  from  ex  open 
var_type  -  char*,  type  of  variable  described 
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“g”  or  “G”  -  global  variables 
“n”  or  “N”  -  nodal  variables 
“e”  or  “E”  -  element  variables 
“m”  or  “M”  -  nodeset  variables 
“s”  or  “S”  -  sideset  variables 

numvars  -  int,  number  of  variables  stored  for  specified  variable  type 
varnames  -  char**,  return  array  of  pointers  to  numvars  variable  names 

Read  Time  Values  for  a  Time  Step 

int  ex  get  time  (exoid,  time  step,  time_value) ; 

exoid  -  int,  returned  value  from  ex  open 
time  step  -  int,  time  step  number  (>=  1) 
timevalue  -  float*,  returned  time  at  the  specified  time  step 

Read  Nodal  Variables 

int  ex_get_nodal_var(exoid,  int  timestep,  nodalvarindex,  numnodes, 
nodal_var_vals); 

exoid  -  int,  returned  value  from  ex  open 

time  step  -  int,  time  step  number  at  which  nodal  variable  values  are  needed  (>=  1) 

nodal  var  index  -  int,  index  of  the  desired  nodal  variable  (>=  1) 

num  nodes  -  int,  number  of  nodes 

nodal  var  vals  -  float*,  returned  array  of  nodal  variables 
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